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Plant breeders have focused on improving plant architec- 
ture as an effective means to increase crop yield. Here, we 
identify the main-effect quantitative trait loci (QTLs) for 
plant shape-related traits in rice (Oryza sativa) and find 
candidate genes by applying whole genome re-sequen- 
cing of two parental cultivars using next-generation se- 
quencing. To identify QTLs influencing plant shape, we 
analyzed six traits: plant height, tiller number, panicle di- 
ameter, panicle length, flag leaf length, and flag leaf width. 
We performed QTL analysis with 178 F 7 recombinant in- 
bred lines (RILs) from a cross of japonica rice line 'SNU- 
SGT and indica rice line 'Milyang23'. Using 131 molecular 
markers, including 28 insertion/deletion markers, we iden- 
tified 11 main- and 16 minor-effect QTLs for the six traits 
with a threshold LOD value > 2.8. Our sequence analysis 
identified fifty-four candidate genes for the main-effect 
QTLs. By further comparison of coding sequences and 
meta-expression profiles between japonica and indica rice 
varieties, we finally chose 15 strong candidate genes for 
the 1 1 main-effect QTLs. Our study shows that the whole- 
genome sequence data substantially enhanced the effi- 
ciency of polymorphic marker development for QTL fine- 
mapping and the identification of possible candidate genes. 
This yields useful genetic resources for breeding high- 
yielding rice cultivars with improved plant architecture. 



INTRODUCTION 

The increasing world population has spawned tremendous 
concerns about impending global food shortages. Nevertheless, 
the per capita increase in global food production has been 
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greater than the increase in population since 1 960, due to enor- 
mous increases in cereal crop production (Trewavas, 2001). 
This dramatic increase in crop yield, known as the Green Revo- 
lution, largely resulted from the development of genetically 
improved high-yielding crop varieties (Hedden, 2003). Among 
the genetic approaches used to improve yield potential, "ideal- 
type breeding" was used to modify plant architecture (morpho- 
logy); for example, the improved rice cultivar IR24 has more 
tillers, reduced plant height and erect leaves. Rice breeders 
have harnessed natural variation in cultivars to modify rice plant 
architecture. For example, incorporating the recessive gene 
sd1 from the Chinese variety Dee-geowoo-gen, reduced plant 
height (Suh and Heu, 1978), and incorporating the photoperiod- 
insensitive gene se1 increased rice adaptability to a wide range 
of latitudes (Sano, 1992; Tamura et al., 1998). 

Although these genetic variants played an important role in 
the marked increase in crop production, the demand for rice 
varieties with higher yield potential and greater yield stability 
has continued to rise due to increasing human population and 
changing global climate (Khush, 2005). To further enhance the 
yield potential of rice over that of existing high-yielding cultivars, 
crop scientists proposed an ideal plant architecture that in- 
cludes fewer tillers with few or no unproductive tillers, more 
grains per panicle, and thick and sturdy stems (Khush, 1995; Li 
et al., 2012). Some of the genes and QTLs regulating architec- 
ture-related traits have been reported. For example, OsSPL14, 
which is regulated by OsmiR156, affects architecture-related 
traits; a point mutation in OsSPL14 perturbs OsmiR156-media- 
ted regulation of OsSPL14 and generates an ideal-type rice 
plant with fewer tillers, increased lodging resistance and en- 
hanced grain yield (Jiao et al., 2010). Marathi et al. (2012) re- 
ported 19 novel QTLs and 15 QTL hotspots affecting plant 
architecture, including yield-related traits, using RILs developed 
from Pusa1266, a semi-dwarf, high-yielding cultivar. Further 
work identified several other QTLs regulating plant shape, in- 
cluding leaf length, leaf width, tiller number and panicle length 
(Farooq et al., 2010; Liu et al., 2010; Wang et al., 2010; 2011). 
Modifying plant architecture for architecture-based breeding will 
require identification of QTLs and their causal genes in different 
genetic backgrounds. 

QTL mapping and cloning using DNA markers are useful 
tools to unlock the genetic basis of complex phenotypic varia- 
tion. Moreover, the application of molecular markers developed 
from QTL analysis enhances breeding efficiency by enabling 
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marker-assisted selection for particular agronomic traits. How- 
ever, traditional approaches for QTL analysis require substan- 
tial time and labor because of the low efficiency of polymorphic 
marker development for high-resolution map construction and 
the difficulties of candidate gene identification. Whole-genome 
sequencing, using next-generation technology can overcome 
these limitations of QTL analysis. A sequencing-based geno- 
typing method uses single nucleotide polymorphisms (SNPs) 
identified by whole-genome sequencing of a mapping popula- 
tion (Huang et al., 2009). This approach could substantially 
reduce the amount of time and effort required for genotyping of 
a mapping population. Recent work also proposed another 
approach, named 'QTL-seq', for rapid identification of QTLs for 
a given phenotype. QTL-seq works by sequencing DNAs from 
two populations, each composed of 20-50 individuals showing 
extreme opposite trait values in segregating progenies (Takagi 
et al., 2013). Although these approaches improve the time- and 
cost-effectiveness of QTL mapping, approaches to increase the 
efficiency of candidate gene identification have not yet been 
proposed. 

In this study, we conducted QTL mapping for six traits related 
to plant architecture in rice and identified 11 main-effect and 16 
minor-effect QTLs. By analysis of whole-genome sequencing 
data and meta-expression profiles, we identified 1 5 strong can- 
didate genes for the 11 main-effect QTLs. The plant-shape 
QTLs and their candidate genes provide useful genetic re- 
sources for breeding of high-yielding rice cultivars. 

MATERIALS AND METHODS 

Plant materials and growth conditions 
To identify QTLs influencing plant shape in rice, we crossed a 
japonica cultivar 'SNU-SG1' with the tongil-type cultivar Mil- 
yang23 with a genetic make-up similar to indica (> 90%) (Yoo 
et al., 2007). Among 338 F 7 RILs developed from the F 2 popu- 
lation by single seed descent, we chose 178 F 7 RILs to conduct 
QTL mapping; these included 40 lines showing extreme pheno- 
typic values for each trait. F 7:8 heterogeneous inbred family- 
near isogenic lines (H-NILs) derived from F 7 RILs were used to 
test the allele effects of candidate gene of plant shape-related 
QTLs (Tuinstra et al., 1997). In 2011, the RIL population and 
two parents were evaluated at the Seoul National University 
Experimental Farm in Suwon (37° N latitude), Korea. The field 
experiment was conducted in a randomized complete block 
design with two replications in two consecutive years (2010 and 
2011) and field management followed normal rice agronomic 
practices. Seeds were sown in the seed bed of the greenhouse 
on April 20 and seedlings were transplanted into the paddy field 
on May 19 at a spacing of 0.3 m by 0.15 m with one seedling 
per hill. Compound fertilizer (21N-17P-17K) was applied in the 
same amount (24 kg/10a) and applied three times: before 
transplanting and at 7 and 10 weeks after transplanting. 

Evaluation of agronomic traits 

For the phenotypic evaluation of the selected 178 F 7 RILs, we 
measured these agronomic traits; plant height, tiller number, 
panicle diameter, panicle length, flag leaf length and flag leaf 
width. All of the traits were evaluated when panicles were fully 
emerged. Panicle diameter (0.5 cm below the panicle base) 
was measured by a digital caliper. Plant height was measured 
from the soil surface to the apex of the main stem. On the main 
stem, flag leaf length was measured from the junction of leaf 
blade and sheath to the leaf apex, and flag leaf width was 
measured at the widest location of the leaf. Tiller number was 



scored when the grains had matured fully. Of all the flowered 
tillers of an individual plant, the longest panicle was measured 
for panicle length. 

Molecular markers for genotyping 

Genomic DNAs were extracted from leaf tissues by a CTAB 
method (Murray and Thompson 1980). Of 131 DNA markers 
used for genotyping, 68 markers are RM-series simple se- 
quence repeat (SSR) markers that were designed by Temnykh 
et al. (2001), and 35 markers are S-series sequence-tagged 
site (STS) markers that were designed based on the sequence 
differences between japonica and indica, using information 
from the Crop Molecular Breeding Lab, Seoul National Univer- 
sity (Chin et al., 2007). The other 28 markers are NGS-based 
insertion/deletion markers designed by sequence comparison 
between the two parent genotypes. To develop these markers, 
we investigated 7 to 1 1 base pairs of unique insertion/deletion 
(InDel) regions by comparing whole-genome sequencing data 
of the two parents, and designed primers to get PCR product (< 
150 bp) containing the InDels. The initial genetic map was con- 
structed with SSR and STS markers, and the next-generation 
based InDel markers were used to further map unlinked re- 
gions in the linkage map. 

Linkage map construction and data analysis 
A molecular linkage map was constructed using Mapmaker 3.0 
(Lander et al., 1987; Lincoln et al., 1993). Distance between mar- 
kers is presented in centiMorgans (cM) using the Kosambi map 
function (Kosambi, 1944), and the order of markers was estab- 
lished using three point linkage analyses. QTL analysis was 
conducted with composite interval mapping implemented in 
Qgene 4.3.10 software (Joehanes and Nelson, 2008). Signifi- 
cance thresholds for QTL detection were determined with 1 000 
permutations, and the threshold of logarithm of odds (LOD) for 
each QTL ranged from 2.24 to 3.03 with 95% confidence. The 
QTLs explaining more than 1 0% of the phenotypic effect were 
defined as main-effect QTLs as described previously (Wang et 
al., 2010). 

High-throughput whole genome re-sequencing 
High-throughput whole-genome sequencing was performed on 
the two parents, SNU-SG1 and Milyang23. The lllumina Ge- 
nome Analyzer II was used to generate short reads and bases 
were called with the Sequence Analysis Software of Pipeline 
version 1.4 of the Genome Analyzer (lllumina, Inc., USA). Ge- 
nomic DNAs were extracted using the QIAquick PCR Purifica- 
tion Kit (Qiagen) from leaf tissues. Full sequencing was pro- 
cessed in the National Instrumentation Center for Environmen- 
tal Management (NICEM) in Seoul National University, Seoul, 
Korea. The sequencing yielded 26X and 31 X fold coverage, as 
well as 9 Gbp and 1 1 Gbp in read length, for SNU-SG1 and 
Milyang23, respectively (Supplementary Tables S1 and S2). 
Then, the 75 bp paired-end reads of the two parents were mapped 
to the japonica cv. Nipponbare reference genome, which was 
used as the reference sequence (International Rice Genome 
Sequencing Program, 2005). Finally, the low-quality bases (Q 
score in scale < 20), and sites with conflicting genotypes among 
reads were excluded, and only the reads aligned to unique 
locations in the reference genome were used for sequence 
construction. The total number of SNPs was 3,202,922, which 
include 163,933 from SNU-SG1 and 3,038,989 from Milyang23. 

Identification of candidate genes for the QTLs 

Of the 27 QTLs, 1 1 main-effect QTLs with high LOD and R 2 
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Table 1 . Correlations of the six agronomic traits analyzed in the F 7 recombinant inbred lines 
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Significance level: ***P < 0.001 , ** < 0.01 , and * < 0.05. 

Ph, plant height; PI, panicle length; Fll, flag leaf length; Flw, flag leaf width; Pd, panicle diameter; Tin, tiller number. 




Fig. 1. Frequency distribution of 
the six agronomic traits in the 
178 F 7 RILs derived from the 
cross of SNU-SG and Milyang23. 
The y-axis of each graph indi- 
cates the number of F 7 RILs. 
Means and ranges of the parents 
are marked at the top of each 
histogram; S, SNU-SG1; M, 
Milyang23. 
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value were used for candidate gene analysis. We identified 
candidate genes harboring mutated sequences such as SNPs 
and InDels by the following process. First, a reference sequen- 
ce (Nipponbare) corresponding to each main-effect QTL was 
transferred from the GRAMENE database (http://www. gra- 
mene.org/) into text in Microsoft Word. The genomic region 
associated with each QTL was selected by confidence interval 
(CI) that was determined by 1 .5-LOD support interval (Dupuis 
and Siegmund, 1999). 1. 5-LOD support intervals give a 95% 
probability of QTL location on the map and were calculated as 
the longest contiguous interval in which the LOD value was 
within 1 .5-fold of its maximum. Coding sequence regions de- 
rived from SNU-SG1 and Milyang23 were interrogated to detect 
mutated sequences using the AutoHotkey macro program 
(http://www.autohotkey.com/). This process has two steps. First, 
the polymorphic sequences (< 15 bp) obtained from whole 
genome sequencing data of the two parents were copied and 
pasted into Microsoft Excel for sequence blast using Auto- 
Hotkey source codes 1 and 2 (Supplementary Fig. S1). Second, 
the sequences were blasted against the Nipponbare genome 
sequence to locate the genomic regions of polymorphic se- 
quences in Microsoft Word using the AutoHotkey source code 
3 (Supplementary Fig. S1). Among the genes harboring mu- 
tated sequences, only the genes for which the mutation caused 
amino acid changes were selected as candidate genes. Arabi- 
dopsis orthologous gene information for the candidate QTL 
genes was obtained from The Arabidopsis Information Re- 



source (TAIR, http://www.arabidopsis.org/). 

Heatmap analysis or microarray data analysis 
We used the information of 1 1 50 Affymetrix array data (platform 
accession number is GPL2025) from the NCBI gene expres- 
sion omnibus (GEO, http://www.ncbi.nlm.nih.gov/ geo/) to ana- 
lyze the mRNA expression pattern of candidate genes in Sup- 
plementary Table S3 (Caoet al., 2012; Jung etal., 2011). From 
the Affymetrix anatomical meta-expression database, we gen- 
erated the log 2 average intensity data of 16 tissues/organs; 
callus, coleoptile, collar, leaf, flag leaf, root, shoot, stem, internode, 
shoot apical meristem (SAM), panicle, anther, stigma/ovary, 
palea/lemma, developing seed, and embryo. Among them, the 
data for stem and palea/lemma were only available for indica, 
the data for coleoptile, collar, internode, stigma/ovary and emb- 
ryo were only for japonica, and the other 9 tissues/organs have 
the data for both. Then, the TIGR Multiple Experiment Viewer 
(MEV) software was used to prepare the Heatmap using the 
anatomical expression data (Supplementary Fig. S2; Table S3) 
(Eisenetal., 1998). 

RNA extraction, reverse transcription (RT) and quantitative 
real-time PCR (qPCR) 

Flag leaf tissues of SNU-SG1 and Milyang23 grown in paddy 
field were homogenized in liquid nitrogen. Total RNA extraction 
and RT-qPCR were performed as previously described (Yoo et 
al., 2009). Briefly, the RT products equivalent to 50 ng of total 
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Table 2. QTL identification for the six agronomic traits by composite interval mapping using 1 78 F 7 recombinant inbred lines 
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a ' b Chromosome number and marker intervals 

"Positive and negative values indicate additive effects contributed by SNU-SG1 and Milyang23 alleles, respectively. 

d Phenotype variation explained by the detected QTL 

e 95% confidence interval for the detected QTL; cM represents centiMorgans. 

f LOD thresholds raised by permutation tests, p = 0.05 



RNA and GoTaq qPCR Master Mix (Promega) were used for 
qPCR reactions using the Light Cycler 480 (Roche) and Roche 
Optical System software was to calculate threshold cycle (DCT) 
value. The rice Ubiquitin gene (Os01 g22490) was used as an 
internal control. The primers used for qPCR are listed in Sup- 
plementary Table S4. 

RESULTS 

Phenotypic variation in F 7 RILs from the Milyang23/SNU- 
SG1 cross 

For QTL analysis of the morphological traits that determine 
plant architecture, we used 178 F 7 RILs from the cross of Mi- 
Iyang23 (indica) and SNU-SG1 (japonica) (Yoo et al., 2007). 
We analyzed six shape-related traits that critically affect lodging 
and yield: plant height, tiller number, panicle diameter, panicle 
length, flag leaf length, and flag leaf width. We measured phe- 
notypic variation in the two parental lines and the RILs (Fig. 1). 



Most of the traits showed approximately normal distribution with 
bidirectional transgressive segregation. We also observed a 
large amount of variation between the two parental cultivars for 
the traits, except for plant height and tiller number. Supplemen- 
tary Table S5 shows descriptive statistics of the six agronomic 
traits. We analyzed the correlation among the traits in the F 7 
RIL population (Table 1) and observed a positive correlation 
between flag leaf length and other traits, except for flag leaf 
width. We also found relatively high positive correlations be- 
tween plant height and panicle length (0.64), and between pa- 
nicle diameter and flag leaf width (0.53). 

Linkage map construction 

A coarse-scale linkage map was initially constructed over the 
whole rice genome using 178 F 7 RILs and 68 SSR and 35 STS 
markers. To construct linkage map for the unlinked regions in 
the initial map, we designed 28 next-generation based InDel 
markers using the whole-genome sequencing data. We used 
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Fig. 2. Chromosomal locations of QTLs de- 
tected by Q-genes in the F 7 RIL population of 
SNU-SG1/Milyang23. Chromosomes are num- 
bered at the top and markers are listed on the 
right of each chromosome. Wave marks in 
chromosome 4 represent the unlinked region. 
Various geometric figures indicate the location 
of peak LODs of the six agronomic traits. 
Chromosome 7 and 12 are not shown be- 
cause we detected no QTLs on these two 
chromosomes. 
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Mapmaker 3.0 to construct the genetic linkage map, which 
consisted of 12 linkage groups and covered 2,451 cM (Fig. 2). 
The average distance between adjacent markers was 14.5 cM, 
which is less than 20 cM, the minimal compulsory level for QTL 
mapping (Lander and Botstein, 1989). 

QTL analysis for the six agronomic traits 
To identify QTLs, we used a LOD threshold value and pheno- 
typic effect (Ft) (see "Materials and Methods"). As a result, we 
identified a total of 27 QTLs for the six traits, with 7.4 to 33.5% 
phenotypic effect (P?) and 2.9 to 15.8 LOD values, across all 
the chromosomes except chromosomes 7 and 12 (Table 2). 
Thus, chromosomes 7 and 12 were excluded from Fig. 2, be- 
cause they did not carry any significant QTL above the thre- 
shold. A QTL with highest LOD (15.8) and phenotypic effect 
(33.5%) was detected for plant height on the long arm of chro- 
mosome 1 (Table 2; Fig. 2). Overall, we detected 3 to 6 QTLs 
for each agronomic trait. Of the 27 QTLs, 1 1 QTLs explained 
more than 1 0% of the phenotypic effect and were defined as 
the main-effect QTLs, which include FII4 for flag leaf length, 
Flw5 and Flw9 for flag leaf width, PI1 and PI9 for panicle length, 



Ph1, Ph3 and Ph9 for plant height, and Tln2, Tln11-1 and 
77n7 1-2 for tiller number (Table 2). Especially, three QTLs (PI9, 
Ph1 and Tln2) explained more than 15% of the phenotypic 
effect, suggesting that these genetic materials provide a valua- 
ble resource for the identification of major genes affecting plant 
architecture. 

QTL hotspots, genomic locations associated with multiple 
traits, may harbor critical regulators of various agronomic traits 
(Neto et al., 2012). Two QTL hotspots were identified on chro- 
mosomes 4 and 9 (grey blocks in Fig. 2). One, in a window of 
around 8 cM on the long arm of chromosome 9, harbors three 
main-effect QTLs (Flw9, Ph9 and PI9 for flag leaf width, panicle 
height, panicle length, respectively) and two minor-effect QTLs 
(Pd9 and F//9for panicle diameter and flag leaf length) (Fig. 2). 
The other, in a window of around 5 cM on chromosome 4, car- 
ries three QTLs, including one main-effect QTL (FII4 for flag leaf 
length) and two minor-effect QTLs (Flw4 for flag leaf width and 
Tln4 tiller number) (Fig. 2). 

To compare the physical locations of the QTLs identified in 
this study and previously reported QTLs, we searched the lite- 
rature for cases in which mapping populations derived from the 
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Table 3. Comparison of the QTLs identified in current and previous studies using the mapping populations derived from crosses of indica and 
japonica cultivars 



QTL of this study QTL of previous studies 





QTL 


Physical 
position (Mb) 


QTL 


Marker region 


Physical 
position (Mb) 


PT a 


PS b 


Mapping parents 
japonica (J)/indica (I) 


References 


Plant 


rnl 


OD.O-Oy.D 


(J / Lptl 7 


OoO-O / 4^ 


oy .0-4^. / 


DM 
DlL 


QQ 

yo 


Nippnbare (J)/Kasalath (I) 


isnirnaru ei al. ^UU4J 


height 






upn \.z 


PfiO££1 QM-1QQ~7 

touoo l-rilvl loo/ 


OQ.U-4U.O 


r 2 


oU I 


iNipponuare ^jyir\io40-ooy \\) 


LIM el al. [£\J 1 1 ) 




rtib 


^o.o-^/.y 


upno.^ 


nIVloo / y-nlvlo4U 


nc o no c 


T 2 


oU l 


NipponDare (jyimo4o-ooy \\) 


Lin et al. \ ) 








Qph6f 


RM162-RM30 


24.0-27.3 


RIL 


226 


z.nongnuiyouo (ixjy 
Xieqingzao B (1) 


Liang etal. (2011) 




Ph9 


16.8-18.1 


Oph9 


RM3912-RM278 


10.8-27.1 


F 2 


301 


Nipponbare (J)/IR1 545-339 (I) 


Linetal. (2011) 


Flag leaf 


FII5 


23.9-27.3 


qFLL-5 


BIN 


27.8-28.6 


RIL 


150 


Nippnbare (J)/931 1 (I) 


Wang etal. (2010) 


length 






QFII5 


RM480-RM334 


27.5-28.6 


RIL 


180 


I RATI 09 (J)/Zhenshan 97(1) 


Yue et al. (2006) 




FII9 


15.9-16.8 


qFII9 


RM24424-RM24434 


16.2-16.5 


F 2 


176 


SN265 (J)/LTH(J) 


Huanq etal. (2010) 


Flag leaf 


Flw4 


28.5-30.7 


qFlw4 


RM17483-RM17486 


31.0-31.2 


RIL F 7 


190 


D50 (J)/HB277(I) 


Chen etal. (2012) 


width 






Qflw4 


RM255-RM349 


28.5-31.6 


RILF 10 


180 


IRAT109 (J)/Zhenshan97(l) 


Yue et al. (2006) 




Flw9 


15.9 


qFLA9 


RM6839B-RM257 


14.5-17.9 


BRIL 


244 


931 1 (l)/Zhenshan97B(l) 


Wang etal. (2012) 


Tiller 


Tln2 


23.0-25.8 


Tln2-2 


RM526-RM425 


26.6-32.2 


SSSL 


78 


IAPAR9 (J)/Hua-Jing-Xian74 (I) 


Liu etal. (2012) 


number 






tp2b 


RM561-RM6318 


19.8-24.4 


BCRIL 


244 


Zhenshan97B(l)/9311(l) 


Liu et al. (2008) 




Tln4 


30.7-33.6 


qTLN4 


RM303-RM349 


27.5-32.4 


RIL-Fn 


254 


Xiushui 79 (J)/C Bao(J) 


Jiang etal. (2012) 








Tin 4-3 


RG143-RG620 


34.2-34.6 


DH d 


123 


Azucena (J)/IR64(I) 


Yanetal. (1998) 




Tln6 


19.4 


pn6 


S6070 


20.1 


F 2 


146 


Junambyeo (J)/IR71 033(1) 


Rahman et al. (2008) 




Tln9 


10.5-12.6 


Tln9 


RM105-RM278 


12.5-19.3 


SSSL 


91 


Hua-Jing-Xian74 (l)/Amol3(l) 


Liu etal. (2012) 


Panicle 


PI1 


36.8-39.6 


qPL-1 


RM472-RM104 


37.8-40.1 


RIL 


187 


IRAT109 (J)/Zheshan97B (I) 


Liu et al. (2008) 


length 


PI4 


32.4 


OPI4b 


RM349-RM280 


32.4-34.9 


BC 2 F 5 


59 


Tarome Molaei (J)/IR64 (I), 


Ahamadi et al. (2008) 



a PT, population type 
b PS, population size 

"Numbers at the end of QTL represent chromosome number 
d DH, double haploid population 



RM341 



RM475 



N2230.1 



RM263 



S20112.6 



RM240 



Line# 



tiller numbers for F 7:8 H-NILs are shown with standard deviation (SD). 




Progeny test 

(F 7B H-NILs) 
Max Min SD 


27 


6 


4.06 


30 


8 


4.37 


21 


2 


4.03 


30 


6 


4.69 


28 


9 


4.20 


22 


5 


3.35 


17 


5 


3.10 



Fig. 3. Development of heteroge- 
neous inbred family-near isogenic 
lines (H-NILs) for TLN2. The rows 
represent the genotypes of seven 
F 7 RILs that contain heterozygous 
genomic regions for the TLN2 QTL. 
Shaded regions indicate chromo- 
somal segments that are homo- 
zygous SNU-SG1 (black), homo- 
zygous Milyang23 (white), and he- 
terozygous Milyang23/H143 (grey). 
40 plants of F 7:8 H-NILs derived 
from each F 7 RIL were investigated 
for progeny tests (on the right). 
Maximum (Max) and minimum (Min) 



crosses of indica and japonica cultivars were used for QTL 
identification. By using flanking DNA markers, we compared the 
physical locations between the QTLs in the current study and 
the previously published QTLs. Among 27 QTLs identified in 
this study, 13 QTLs were in accord with the chromosomal re- 
gions carrying the published QTLs (Table 3). The remaining 14 
QTLs did not overlap with the published QTLs, indicating that 



these are potentially novel QTLs for the traits. Notably, the 
QTLs for panicle diameter are all newly discovered. Among the 
new 14 QTLs, 6 QTLs are the main-effect QTLs with above 
10% of phenotypic effect, including FII4, Flw5, PI9, Ph3, Tln11- 
1 and 77nn-2(Table2). 
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Table 4. List of 1 1 strong candidate genes for the main-effect QTLs obtained by sequence analysis of whole genome re-sequencing data for 
the two parent cultivars, Milyang23 {indica) and SNU-SG1 (japonica) 



QTL Type 3 Gene ID Mutation" Description Ortholog 0 



PI1, Ph1 


Milyang23 


Os01g68510 


133/601 


Peptide transporter PTR2 


AT1G69870 (NRT1.7) 


Tln2 


SNU-SG1 


Os02g40100 


398/791 


DUF869 domain containing protein 


AT2G23360 


Tln2 


SNU-SG1 


Os02g40190 


33/240 


Receptor kinase 




Tln2 


SNU-SG1 


Os02g40200 


485/1033 


Receptor-like protein kinase precursor 




Tln2 


SNU-SG1 


Os02g40240 


675/1002 


Receptor-like protein kinase precursor 




Tln2 


SNU-SG1 


Os02g40280 


273/370 


Piwi domain containing protein 


AT1G69440(AGO7) 


Tln2 


SNU-SG1 


Os02g41610 


642/745 


Expressed protein, Coiled coil 




FII4 


SNU-SG1 


Os04g51450 


35/331 


Glycosyl hydrolases family 16 


AT2G18800 (XTH21) 


FII4 


SNU-SG1 


Os04g53120 


6A.S. 


XA1 , NB-ARC domain containing protein 




Flw5 


SNU-SG1 


Os05g31230 


5A.S. 


N-acetyl transferase ESC01 


AT4G31400 


Tln11-1 


Milyang23 


Os11g18194 


560/761 


Cycloartenol synthase 


AT2G07050 (CAS1) 



a Parent type harboring mutant alleles, which were determined by sequence comparison with reference genome japonica cv. Nipponbare. 

"The position of amino acid in which premature stop codon occurred/total protein length, A.S. amino acid substitutions. 

■Accession number of Arabidopsis orthologs were obtained from The Arabidopsis Information Resource (TAIR, http://www.arabidopsis.org/). 



Table 5. Seven strong candidate genes for the main-effect QTLs obtained by expression profile analysis of Affymetrix array data 



Gene ID a 


QTL Name" 


Description 




Log 2 fold change (japonica/ indicaf 




Flag leaf 


Root 


Shoot 


SAM d 


Panicle 


Anther 


Os09g26530 


Flw9, Ph9, PI9 


Ser/Cys_Pept_Trypsin-like 


2.9 


2.3 


1.6 


2.6 


2.4 


2.1 


Os09g26554 


Flw9, Ph9, PI9 


Ser/Cys_Pept_Trypsin-like 


3.4 


2.5 


2.0 


3.0 


3.2 


3.5 


Os03g26080 


Ph3 


Nucleoside triphosphatase 


2.3 


1.1 


1.0 


1.3 


1.1 


1.0 


Os11g20239 


Tln11-1 


Hydroxyproline-rich glycoprotein 
DZ-HRGP precursor 


2.3 


2.0 


2.3 


2.8 


1.9 


2.0 


Os02g40190 


Tln2 


Receptor kinase 


2.4 


1.1 


1.9 


0.8 


0.7 


1.1 


Os02g40200 


Tln2 


Receptor-like protein kinase precursor 


2.4 


1.4 


2.3 


2.1 


0.9 


0.8 


Os02g40240 


Tln2 


Receptor-like protein kinase precursor 


4.4 


1.9 


2.9 


1.0 


1.5 


1.5 



a Locus identifier (ID) annotated by Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/) team 

b QTL name queried from Gramene QTL database (www.gramene.org/qtl/) 

c Log 2 fold change of japonica/ indica in Affymetrix anatomical meta-expression profile 

d SAM represents shoot apical meristem. 



Sequencing analysis for candidate gene identification 
We carried out whole-genome sequencing of the two parenta- 
cultivars, SNU-SG1 and Milyang23, to obtain a cultivar-specific 
nucleotide polymorphism data set (see "Materials and Me- 
thods"). The polymorphic nucleotide sequences, including In- 
Dels, were obtained by comparing the whole genome se- 
quences of SNU-SG1 and Milyang23 with the reference Nip- 
ponbare sequence using Pipeline version 1.4 (see "Materials 
and Methods"). Of the total 363,845 SNPs identified by compar- 
ing with the reference sequence, 330,317 were identified from 
Milyang23, and 33,528 from SNU-SG1. Sequence variation 
between indica- and yapon/ca-type subspecies contributed to 
the large number of SNPs in Milyang23. For identification of 
candidate genes associated with the QTLs, we analyzed all the 
SNPs (including InDels) derived from SNU-SG1, but only In- 
Dels capable of impairing protein function were considered for 
Milyang23 to filter out non-causal subspecies sequence varia- 
tion between indica and japonica genomes. 

To perform candidate gene analysis for 1 1 main-effect QTLs, 
we analyzed the genomic region that was selected based on 



confidence interval for each QTL peak (Table 2). Confidence 
intervals were determined by 1 .5-LOD support interval, which 
provides a QTL coverage probability of approximately 95% 
(Dupuis and Siegmund, 1999). The confidence intervals of the 
QTLs detected in the integrated map using composite interval 
mapping are given in Table 2. We determined the specific loca- 
tions of polymorphic sequences in the individual genes by 
BLAST to the Nipponbare reference sequence. By these cod- 
ing sequence analyses, we identified 54 genes with 656 nuc- 
leotide sequence polymorphisms in their open reading frames 
within the 1 1 main-effect QTL regions (data not shown). Among 
them, only the 308 SNPs that caused amino acid variation were 
considered as candidate genes. Thus, by analysis of the whole- 
genome sequencing data, we isolated 54 candidate genes 
carrying at least one amino acid polymorphism for the six traits 
(Supplementary Table S3). 

Strong candidate genes for the main-effect QTLs 

Among the 54 candidate genes for 1 1 main-effect QTLs, we 
further isolated strong candidate genes by analyzing their ami- 
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Fig. 4. Comparison of efficiencies of poly- 
morphic marker development between next- 
generation sequencing based InDel (NID) 
markers and conventional SSR and STS 
markers. The numbers in parentheses 
represent the numbers of polymorphic or 
non-polymorphic markers out of total DNA 
markers that were tested for marker devel- 
opment. Polymorphic markers represent 
the markers obviously distinguishable in 
PCR product size by 3% agarose gel elec- 
trophoresis. NID markers were used to 
construct a linkage map for unlinked re- 
gions in the initial map. 
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Fig. 5. Genomic locations of two main- 
effect QTLs and their candidate genes. (A, 
B) Os01g68510 (A) and Os11g18194 (B) 
genes located at two main-effect QTLs, 
Ph1 and Tln11-1, respectively. (C, D) Gene 
structures and natural variations of 
Os01g68510 (C) and Os11g18194 (D). 
Sequence variations between SNU-SG1 
and Milyang23 are indicated. Black dots 
represent the absence of the correspond- 
ing bases. 



M23 
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no acid sequences and meta-expression profiles. First, by pro- 
tein coding sequence analysis, we selected 1 1 genes as strong 
candidates, because they harbor potential causal SNPs leading 
to a premature stop codon or frameshift in their amino acid 
sequences (Table 4). These 11 strong candidate genes inclu- 
ded candidates for all agronomic traits examined, except for 
panicle diameter, with one to seven genes for each trait. These 
included single genes for panicle length, plant height and flag 
leaf width, two genes for flag leaf length, and seven genes for 
tiller number (Table 4). One candidate gene (Os01g68510), 
encoding peptide transporter 2 (PTR2), was identified for two 
QTLs related to panicle length and plant height on the distal 
part of chromosome 1 (Table 4). The largest effect QTL Ph1 for 
plant height was also detected with 33.5% P? value in this loca- 
tion on chromosome 1 (Tables 2 and 3). Meanwhile, QTL effect 
can be also caused by differences in gene expression between 



alleles due to mutations in regulatory regions including promo- 
ters and introns (Meuwissen, 2010). Thus, we next analyzed 
meta-expression profiles of the 54 candidate genes using Af- 
fymetrix array data (Supplementary Fig. S2), which provide the 
indica or yapon/ca-specific gene expression (Cao et al., 2012; 
Jung et al., 201 1). As we used the indica (Milyang23) and japo- 
nica (SNU-SG1) cultivars to create the QTL mapping popula- 
tion, the genes expressed differentially in the indica and japoni- 
ca rice cultivars would be considered as strong candidates for 
the corresponding traits. By Heatmap analysis, we found that 
seven genes are differentially expressed in at least one tissue 
type with more than two-fold differences between indica and 
japonica (Table 5). Interestingly, all the probes corresponding to 
the seven genes were expressed preferentially in the flag 
leaves of japonica rice. To test whether these genes are diffe- 
rentially expressed in SNU-SG1 and Milyang23, we performed 
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RT-qPCR analysis using flag leaves. Of the seven genes, five 
genes were preferentially expressed in SNU-SG1 with more 
than two-fold difference compared to Milyang23 (Supplementa- 
ry Fig. S3), in accordance with meta-expression profile. How- 
ever, we failed to amplify the other two genes (Os02g40190 
and Os02g40240), possibly due to very low transcript abun- 
dance. This result strongly supports the idea that at least five 
genes identified by meta-expression profile analysis are strong 
candidate genes for the corresponding QTLs. Notably, three 
genes (Os02g40190, Os02g40200 and Os02g40240) for Tln2 
were identified as strong candidates by both coding sequence 
and meta-expression analyses, and encode known receptor- 
like protein kinase precursors (Tables 4 and 5). Thus, taken 
together our data provide a total of 1 5 strong candidate genes 
by two different analyses. 

We further investigated the orthologs of these candidate 
genes in Arabidopsis thaliana, using Gramene (http://www. 
gramene.org/) and TAIR (http://www.arabidopsis.org/) data- 
bases. We found six Arabidopsis orthologs in which four genes 
have been functionally characterized (Table 4). In PI1/Ph1 QTL, 
At1g69870 encodes a low affinity nitrate transporter (NRT1.7), 
which functions in nitrate remobilization in leaf (Fan et al., 2009). 
In Tln2 QTL, At1g69440 encodes ARGONAUTE 7, a member 
of the ARGONAUTE family, which functions in the regulation of 
developmental timing (Adenot et al., 2006). In FII4 QTL, At2g 
18800 encodes Xyloglucan endotransglycosylase/hydrolases 
(XTHs) which functions in carbohydrate metabolism (Paul and 
Janet, 1999). In Tln1-1 QTL, At2g07050 encodes cycloartenol 
synthase 11 (Corey et al., 1993) which functions in the biosyn- 
thesis of brassinosteroids (Asami and Yoshida, 1999). This 
identification of promising candidate genes provides useful 
genetic resources for the further isolation of the causal gene 
associated with each QTL. 

Validation of QTL regions using H-NILs 
To further validate the QTL regions associated with 15 strong 
candidate genes, we developed H-NILs that segregate for the 
phenotypes of the corresponding traits (Tuinstra et al., 1997). 
The genotypic and phenotypic analyses of the progeny H-NILs 
allow identification of the candidate genes underlying QTLs 
(Routaboul et al., 2012). In 178 F 7 RILs, we selected heteroge- 
neous lines for at least one of the flanking marker loci in the 
QTL regions (Fig. 3). By phenotypic analysis using 40 F 7:8 indi- 
vidual plants derived from each heterogeneous line, we found 
several H-NILs for six main-effect QTLs that harbor strong can- 
didate genes (Supplementary Table S6). Most of the H-NILs, 
including F 7 RIL #72-derived ones which showed highest varia- 
tion in tiller number, displayed transgressive segregation (Fig. 3, 
Supplementary S4). These support that the target genes are 
localized within the corresponding QTL regions, enhancing the 
fidelity of strong candidate genes. Further linkage analysis of 
the markers and traits in the H-NIL mapping populations may 
lead to isolation of the genes underlying the QTLs. 

DISCUSSION 

Next-generation sequencing data enhance the efficiency of 
QTL analysis 

Agronomically important crop traits are generally controlled by 
multiple genes (QTLs) that produce continuous variation in 
phenotype. Conventional QTL analysis remains time-consu- 
ming and labor-intensive, mainly because it requires develop- 
ment of polymorphic markers for linkage analysis. To overcome 
this limitation, several research groups have recently introduced 



next-generations sequencing technology as a new strategy for 
genetic mapping. For example, high-throughput genotyping of 
recombinant populations can be performed by whole-genome 
sequencing of individual RILs (Huang et al., 2009). This se- 
quence-based genotyping has proven faster for genetic map 
construction and more accurate for determination of recombina- 
tion breakpoints than marker-based genotyping methods (Wang 
et al., 2011). Another group performed whole genome re- 
sequencing of 50 rice cultivars and investigated genome-wide 
variation in high-quality candidate SNPs obtained from re- 
sequencing data for each cultivar (Xu et al., 2012). In this study, 
we modified conventional mapping methods by whole-genome 
sequencing of the two parental cultivars to improve QTL map- 
ping efficiency. To develop polymorphic DNA markers for lin- 
kage analysis, we generated InDel markers using whole- 
genome sequencing data; these markers were designed based 
on nucleotide sequence polymorphisms between the two par- 
ents. The InDel markers showed much higher efficiency 
(74.5%) in polymorphic marker development compared to con- 
ventional SSR (18.1%) and STS (21.1%) markers, which were 
designed from sequence polymorphisms between japonica and 
indica genomes (Fig. 4). These results indicate that application 
of a whole-genome sequencing approach to the two parental 
cultivars substantially enhanced the efficiency of QTL mapping 
with relatively lower cost for whole-genome sequencing com- 
pared to high-throughput approaches. 

High-resolution mapping to identify candidate genes for the 
detected QTLs is another time-consuming and laborious step in 
conventional QTL analysis. Thus, instead of performing fine 
mapping to narrow down the QTL region to less than 100 kb, 
we directly analyzed the genomic regions for the QTLs that 
were selected based on a 1 .5-LOD confidence interval (Table 
2), and identified candidate genes for each QTL of interest. Our 
sequence analysis found an average of five candidate genes 
harboring 23 causal SNPs in the confidence interval-based 
QTL regions by comparing genomic sequences of SNU-SG1 
and Milyang23 with the reference Nipponbare. In conventional 
QTL analysis, fine mapping has to be conducted to narrow 
down the QTL region to 50 to 1 00 kb through development of 
large numbers of polymorphic DNA markers in each QTL re- 
gion as well as development of mapping populations such as 
heterogeneous inbred family-near isogenic lines (Tisne et al., 
2008) and backcross inbred lines (Li et al., 2011). Whole- 
genome sequencing-based QTL analysis enabled us to isolate 
candidate genes for all the QTLs identified from initial mapping 
population at the same time, while conventional QTL ap- 
proaches need to develop a high resolution-mapping popula- 
tion for each QTL to find strong candidate gene(s). Furthermore, 
the combination of whole-genome sequencing-based QTL and 
meta-expression profile analyses enable us to isolate candidate 
genes for the previously reported QTLs by sequencing of the 
two parents that were used to construct mapping populations. 
Therefore, we propose that this advanced QTL analysis ap- 
proach, sequencing two parent cultivars, is an efficient method 
for candidate gene identification, as it skips the fine-mapping 
procedure. 

Comparison of the identified QTLs with previously 
reported QTLs for agronomic traits 

Some QTLs were identified here and in previous studies of 
plant architecture that examined traits such as tiller number and 
panicle diameter. For example, the four QTLs, Tln2, Tln4, Tln6 
and Tln9, which influence tiller number, shared a similar ge- 
nomic location with QTLs identified in previous studies (Jiang et 
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al., 2012; Liu et al., 2012; Qu et al., 2012; Rahman et al., 2008; 
Yan et al., 1998) (Table 3). This consistency in the QTL regions 
associated with the tiller number among independent studies 
indicates that the locations of alleles for the trait are conserved 
across different genetic and environmental backgrounds. By 
contrast, the other two QTLs, Tln11-1 and Tln11-2, have not 
been previously identified. All the QTLs we identified for panicle 
diameter (Pd3, Pd6, Pd8, Pd9 and Pd1 1) are novel, with no 
correspondence to previously reported QTLs for this trait (Table 
2), indicating that these QTLs may be a potentially new set of 
alleles specific to this genetic background. 

Analysing candidate genes and inferring their functions 
from Arabidopsis homologs 

By the whole-genome sequencing data analysis, we found 54 
candidate genes in all the detected QTL regions related to six 
agronomic traits (Supplementary Table S3). By further se- 
quence and meta-expression profile analyses, We designated 
15 out of 54 genes as strong candidates for 1 1 QTLs. Of the 15 
strong candidate genes, 11 genes identified by amino acid 
sequence analysis carry polymorphic SNPs that lead to a fra- 
meshift or premature stop codon in the coding sequence (Table 
4). For example, in the PI1 and Ph1 QTL region (Fig. 5A), five 
genes containing polymorphic SNPs in their coding sequences 
were initially found by sequence analysis (Supplementary Ta- 
ble S3). Of these SNPs, a 1-bp deletion in the coding region of 
PTR2 occurred in the Milyang23 allele, leading to a premature 
stop codon (Fig. 5C), whereas simple SNPs, which do not 
cause premature stop codon, were found for the other four 
genes. These suggest that PTR2 is the strongest candidate 
gene for the PI1 and Ph1 QTL (Supplementary Table S3). In 
the same way, we identified 1 0 strong candidate genes for the 
other five agronomic traits (Figs. 4B and 4D; Supplementary 
Table S3). 

In addition to mutations in coding sequences, nucleotide mu- 
tations in the promoter, 5'-UTR, 3'-UTR and intron of many 
eukaryotic genes affect morphological traits by regulating gene 
expression at transcriptional or post-transcriptional levels (Lu et 
al., 2008; Rose et al., 2008; Yu et al., 2007). For example, 
TAC1, which affects tiller angle, is highly expressed in the Mi- 
ca rice line 'IR24' which displays a relatively spread-out plant 
architecture, but the expression of TAC1 is low in the intro- 
gressed line, IL55, derived from the japonica rice 'Asominori' 
which displays extremely erect tillers (Yu et al., 2007). This 
differential expression results from a nucleotide change in an 
intron, which affects 3'-UTR splicing. Recently, Jung et al. (2013) 
analyzed publicly available Affymetrix array data and identified 
490 genes preferentially expressed in japonica rice and 104 
genes in indica rice. Some of the genes that are preferentially 
expressed in each subspecies involve development of morpho- 
logical traits such as plant height, panicle, culm and root (Jung 
et al., 2013). Thus, using the Affymetrix array data, we ana- 
lyzed the expression patterns of 54 candidate genes in the 
indica and japonica rice samples of various tissue/organs and 
found seven genes that are differentially expressed between 
japonica and indica; this /'apon/ca-preferred expression was 
confirmed in SNU-SG1 by RT-qPCR analysis (Supplementary 
Fig. S3). Notably, three of the differentially expressed genes 
(Os02g40190, Os02g40200 and Os02g40240) were identified 
by both the sequence and expression analyses (Tables 5 and 
Supplementary Fig. S3). Together with 1 1 candidate genes that 
were identified by sequence analysis, an additional four genes 
from the meta-expression analysis between japonica and indica 
cultivars would be strong candidates for the causal genes for 



the corresponding QTLs. 

The functions of the uncharacterized genes in rice can gen- 
erally be predicted from the functions of orthologous genes in 
Arabidopsis, as proteins in these species generally have con- 
served molecular and biochemical functions. For example, the 
functions of CO and FT, flowering promoters in photoperiod 
pathway, are well conserved between Arabidopsis and rice and 
most likely in other species as well (Hayama and Coupland, 
2004). Among the 15 strong candidate genes proposed in this 
study, the functions of four Arabidopsis orthologs have been 
reported (Corey et al., 1993; Fan et al., 2009; Maris et al., 2009) 
(Table 4). The Arabidopsis nitrate transporter NRT1.7 (Fan et 
al., 2009), an ortholog of PTR2, plays an important role in 
source-to-sink remobilization of nitrate. In addition, a rice gene, 
Short Panicle 1 (SP1) encoding a putative peptide transporter 
(PTR) has been reported to regulate inflorescence branch elon- 
gation; the knock-down sp1 mutant is defective in rice panicle 
elongation, producing short-panicle phenotypes (Li et al., 2009). 
These also indicate that PTR2 is a promising candidate gene 
related to the two QTLs detected for panicle length and plant 
height in this study. Another candidate gene Os11g18194, 
located in the region of Tln11-1 (Fig. 5B), is a rice ortholog of 
Arabidopsis CAS1, encoding a cycloartenol synthase, which is 
involved in the biosynthesis of brassinosteroids (Corey et al., 
1993). Mutations in the genes involved in the BR biosynthesis 
pathway modify tiller number (Fujii et al., 1991), leaf size (Mori- 
naka et al., 2006) and yield per plant (Choe et al., 2001). More- 
over, increased levels of brassinosteroids cause a significant 
effect, increasing tiller number and seed weight (Wu et al., 
2008). These suggest that OsCASI is a promising candidate 
gene for Tln11-1 related to tiller number. Lastly, another candi- 
date gene Os04g51450, located in the region of FII4, has an 
orthologous gene At2g18800 (XTH21), encoding a member of 
the xyloglucan endotransglucosylase/ hydrolases (XTH), which 
is involved in cell loosing process at the preliminary stage of cell 
elongation in shoot and root meristems (Liu et al., 2007). 
Genes related to cell elongation, such as ANGUSTIFOLIA in 
Arabidopsis, play a role in leaf longitudinal development (To- 
mohiko et al., 1996). This supports the possible involvement of 
Os04g51450 in FII4 QTL regulating flag leaf length. 

In this study, we developed an improved approach for DNA 
marker development and candidate gene identification by whole- 
genome sequencing of the two parent cultivars of our RILs. 
This approach substantially increased the efficiency of linkage 
analysis by enabling marker development by simple compari- 
son of sequence data from the two parents. In addition, the 
application of whole-genome sequencing to candidate gene 
analysis reduced the cost and time for candidate gene identifi- 
cation by skipping the high-resolution mapping step. By QTL 
analysis, we detected a total of 27 QTLs for six agronomic traits. 
The 15 novel QTLs, two QTL hotspots and 15 strong candidate 
genes identified in this study provide valuable resources for 
molecular breeding of new high-yielding plant varieties in rice. 

Note: Supplementary materials are available on the Molecules 
and Cells website (www.molcells.org). 
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